1 Gaussian Sequence Model

Recall Gaussian sequence model ${\displaystyle X\sim {\mathcal{N}}_d(\theta ,I_d)}$. The goal is to estimate ${\displaystyle \theta \in {\mathbb{R}}^d}$ via ${\displaystyle \delta (X)}$ with low MSE: $$\mathrm{MSE}(\theta ;\delta )={\mathbb{E}}_{\theta }||\theta -\delta (X)|{|}^2.$$
The model is more general than it appears. Like if ${\displaystyle X_1,\cdots ,X_n\stackrel{\mathrm{i}.\mathrm{i}.\mathrm{d}}{\sim }{\mathcal{N}}_d(\theta ,{\sigma }^2{I}_d)}$ for known ${\displaystyle {\sigma }^2>0}$, we could make a sufficiency reduction to obtain $$Z=\frac{1}{\sigma \sqrt{n}}\sum _{i=1}^n{X}_i\sim {\mathcal{N}}_d(\theta ,I_d).$$

1.1 Bayes Estimators

If we introduce Bayesian prior ${\displaystyle {\theta }_i\stackrel{\mathrm{i}.\mathrm{i}.\mathrm{d}}{\sim }\mathcal{N}(0,{\tau }^2)}$, then Bayes estimator is ${\displaystyle \frac{{\tau }^2}{1+{\tau }^2}X}$.
We can think of this as a tuning parameter for a generic linear shrinkage estimator $${\delta }_{\zeta }(X)=(1-\zeta )X,$$ where ${\displaystyle \zeta \in [0,1]}$ is in effect a tuning parameter we will call shrinkage parameter. Taking ${\displaystyle \zeta =(1+{\tau }^2{)}^{-1}}$ corresponds to the Bayes estimator.
If we aren't sure which ${\displaystyle \zeta }$ to use, (like have a priori uncertainty about ${\displaystyle {\tau }^2}$), we can use hierarchical Bayes: $$\delta (X)=(1-\mathbb{E}[\zeta |X])X={\delta }_{{\hat{\zeta }}_{\mathrm{Bayes}(X)}}(X),$$ so we are in effect estimating ${\displaystyle \zeta }$ from the whole data set and the plugging it in as a data-adaptive tuning parameter.

If ${\displaystyle d\ge 3}$ then the UMVUE for ${\displaystyle \zeta }$ is $${\hat{\zeta }}_{\mathrm{UMVU}}(X)=\frac{d-2}{||X|{|}^2},$$ based on the fact that $$Y\sim {\chi }_d^2=\mathrm{Gamma}(\frac{d}{2},2),d>2\Rightarrow \mathbb{E}\left[\frac{1}{Y}\right]=\frac{1}{d-2}.$$
Plugging in ${\displaystyle {\hat{\zeta }}_{\mathrm{UMVU}}}$ results in an estimator called James-Stein estimator: $${\delta }_{\mathrm{JS}}(X)=(1-\frac{d-2}{||X|{|}^2})X={\delta }_{{\hat{\zeta }}_{\mathrm{UMVU}}}(X).$$

1.2 James-Stein Paradox

While James-Stein estimator can be motivated as an empirical Bayes estimator, it is good even without making any Bayesian assumptions at all.

For ${\displaystyle d\ge 3}$, the estimator ${\displaystyle X}$ is actually inadmissible as an estimator of ${\displaystyle \theta }$: $$\mathrm{MSE}(\theta ,{\delta }_{\mathrm{JS}})<\mathrm{MSE}(\theta ,X),\mathrm{\forall }\theta \in {\mathbb{R}}^d.$$

This is surprising because it not only beats the UMVUE on average, but every fixed value of ${\displaystyle \theta }$.

In fact, if we use an estimator shrinking towards ${\displaystyle {\theta }_0\in {\mathbb{R}}^d}$: ${\displaystyle \tilde{\delta }(X)={\theta }_0+(1-\frac{d-2}{||X-{\theta }_0|{|}^2})(X-{\theta }_0)}$. This also dominates ${\displaystyle {\delta }_0}$ because it is just the James-Stein estimator we'd get if we substitute $$Y=X-{\theta }_0\sim {\mathcal{N}}_d(\mu ,I_d),\mu =\theta -{\theta }_0.$$
By the translation invariance of Gaussian location model, James-Stein estimator for ${\displaystyle \mu }$ using ${\displaystyle Y}$ also dominates ${\displaystyle {\hat{\mu }}_0(Y)=Y}$, i.e. ${\displaystyle {\delta }_0(X)={\hat{\mu }}_0+{\theta }_0=X}$ for ${\displaystyle \theta }$.

1.3 Linear Shrinkage Estimators

Even without introducing a Bayesian prior for ${\displaystyle \theta }$, we can motivate our linear shrinkage estimator purely from the perspective of trading bias for a reduction in variance. Calculate MSE: $$\begin{array}{rl}{\mathbb{E}}_{\theta }[(\theta -{\delta }_i(X){)}^2]& =({\theta }_i-{\mathbb{E}}_{\theta }(1-\zeta )X_i{)}^2+{\mathrm{Var}}_{\theta }(1-\zeta )X_i\\ & =(\zeta {\theta }_i{)}^2+(1-\zeta {)}^2,\\ \Rightarrow \mathrm{MSE}(\theta ;\delta )& ={\zeta }^2||\theta |{|}^2+d(1-\zeta {)}^2.\end{array}$$
So let $$\begin{array}{rl}& 0=\frac{\mathrm{d}}{\mathrm{d}\zeta }\mathrm{MSE}(\theta ;\delta )=2\zeta ||\theta |{|}^2-2(1-\zeta )d\\ \Rightarrow & {\zeta }^{\ast }(\theta )=\frac{d}{d+||\theta |{|}^2}.\end{array}$$
This looks similar to ${\displaystyle \frac{1}{1+{\tau }^2}}$ (Bayes-optimal ${\displaystyle \zeta }$ under the Gaussian prior)

2 SURE

Now consider the multivariate version. For a function ${\displaystyle h:{\mathbb{R}}^d\to {\mathbb{R}}^d}$, define Jacobian matrix ${\displaystyle Dh\in {\mathbb{R}}^{d\times d}}$: $$(Dh(x){)}_{ij}=\frac{\partial h_i}{\partial x_j}(x)$$ and the Frobenius norm ${\displaystyle A\in {\mathbb{R}}^{d\times d}}$ as $$||A|{|}_F=\sqrt{\sum _{i,j}{A}_{ij}^2}.$$

Apply Stein's lemma to ${\displaystyle h(x)=x-\delta (x)}$, then $$\begin{array}{rl}\mathrm{MSE}(\theta ;\delta )& ={\mathbb{E}}_{\theta }||\delta (X)-\theta |{|}^2={\mathbb{E}}_{\theta }||X-h(X)-\theta |{|}^2\\ & ={\mathbb{E}}_{\theta }||X-\theta |{|}^2+{\mathbb{E}}_{\theta }||h(X)|{|}^2-2{\mathbb{E}}_{\theta }[(X-\theta {)}^{\mathrm{T}}h(X)]\\ & ={\sigma }^{2}d+{\mathbb{E}}_{\theta }||h(X)|{|}^2-2{\mathbb{E}}_{\theta }\mathrm{tr}(Dh(X)).\end{array}$$
(Note that we assume ${\displaystyle X\sim {\mathcal{N}}_d(\theta ,{\sigma }^2{I}_d)}$ now.)
So if ${\displaystyle {\sigma }^2}$ is known, we obtain the unbiased estimator $$\begin{array}{}\text{(2.3)}& \widehat{\mathrm{MSE}}(X)={\sigma }^{2}d+||h(X)|{|}^2-2{\sigma }^2\mathrm{tr}(Dh(X)).\end{array}$$
We call it Stein Unbiased Risk Estimator (SURE).

3 Risk of the James-Stein Estimator

Now we calculate the risk of James-Stein estimator ${\displaystyle {\delta }_{\mathrm{JS}}(X)=(1-\frac{d-2}{||X|{|}^2})X}$.

• First assume ${\displaystyle {\sigma }^2=1}$. Like above, denote $$h(X)=\frac{d-2}{||X|{|}^2}X\Rightarrow ||h(X)|{|}^2=\frac{(d-2{)}^2}{||X|{|}^2},$$ and by ${\displaystyle h(X)}$ we have $$Dh(X{)}_{ii}=(d-2)\frac{\partial }{\partial X_i}\frac{X_i}{\sum _{j=1}^d{X}_j^2}=(d-2)\frac{||X|{|}^2-2X_i^2}{||X|{|}^4},$$ then $$\mathrm{tr}(Dh(X))=(d-2)\frac{d||X|{|}^2-2||X|{|}^2}{||X|{|}^4}=\frac{(d-2{)}^2}{||X|{|}^2}.$$
  Plug in (2.3), $$\widehat{\mathrm{MSE}}(X)=d+\frac{(d-2{)}^2}{||X|{|}^2}-2\frac{(d-2{)}^2}{||X|{|}^2}=d-\frac{(d-2{)}^2}{||X|{|}^2}.$$
  Taking expectations: $$\mathrm{MSE}(0;{\delta }_{JS})=d-(d-2{)}^2\frac{1}{d-2}=2.$$
  (Fact mentioned here) The total MSE does not rise with ${\displaystyle d}$!
  On the other hand, suppose ${\displaystyle ||\theta |{|}^2\to \mathrm{\infty }}$, then ${\displaystyle {\mathbb{E}}_{\theta }||X|{|}^2\to \mathrm{\infty }}$, so ${\displaystyle \frac{(d-2{)}^2}{{\mathbb{E}}_{\theta }||X|{|}^2}}$ will be driven to ${\displaystyle 0}$.

• For general ${\displaystyle {\sigma }^2}$, similarly, James-Stein estimator is ${\displaystyle (1-\frac{(d-2){\sigma }^2}{||X|{|}^2})X}$, then ${\displaystyle \widehat{\mathrm{MSE}}(X)={\sigma }^{2}d-{\sigma }^4\frac{(d-2{)}^2}{||X|{|}^2}}$. Then ${\displaystyle \mathrm{MSE}(0;{\delta }_{\mathrm{JS}})=2{\sigma }^2}$.